Mon. Not. R. Astron. Soc. 000, 1-13 (2013) Printed 18 March 2013 (MN IMbX style file v2.2) 



Three-dimensional Keplerian orbit-superposition models of the 
nucleus of M31 



o 

<N 



C. K. Brown and J. Magorrian* 

Rudolf Peierls Centre for Theoretical Astrophysics, 1 Keble Road, Oxford 0X1 3NP 



18 March 2013 



o 

u 

6 



> 

in 

cn 
\o 
cn 

cn 
o 
cn 



ABSTRACT 

We present three-dimensional eccentric disc models of the nucleus of M3 1 , modelling the 
disc as a linear combination of thick rings of massless stars orbiting in the potential of a 
central black hole. Our models are nonparametric generalisations of the parametric models of 
Peiris & Tremaine. The models reproduce well the observed WFPC2 photometry, the detailed 
line-of-sight velocity distributions from STIS observations along PI and P2, together with the 
qualitative features of the OASIS kinematic maps. 

We confirm Peiris & Tremaine's finding that nuclear discs aligned with the larger disc 
of M31 are strongly ruled out. Our optimal model is inclined at 57° with respect to the line 
of sight of M31 and has position angle P.A.= 0/ + 90° = 55°. It has a central black hole of 
mass M, 1.0 x 10^ M©, and, when viewed in three dimensions, shows a clear enhancement 
in the density of stars around the black hole. The distribution of orbit eccentricities in our 
models is similar to Peiris & Tremaine's model, but we find significantly different inclination 
distributions, which might provide valuable clues to the origin of the disc. 
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1 INTRODUCTION 

As the nearest spiral galaxy to our own, M31 allows the study of 
a galactic centre in unmatched detail. A particularly striking fea- 
ture of high-resolution V- and /-band photometry of the central 
few arcseconds of M31 is its double nucleus (Light et al. 1974; 
Lauer et al. 1993; Bacon et al. 1994; King et al. 1995): there are 
two peaks in surface brightness, the brighter of which (known as 
PI) is extended and offset ~ Of'5 from the fainter peak (known as 
P2), which is elongated and close to the photometric centre of the 
bulge. Early spectroscopic observations resolved steep gradients in 
the stellar rotation velocities and a prominent peak in the stellar 
velocity dispersion, which hints at the presence of a massive black 
hole (Dressier & Richstone 1988; Kormendy 1988; Bacon et al. 
1994; van der Marel et al. 1994). 

Tremaine (1995) put forward an elegant explanation for these 
observations: the nucleus is a massive disc of red stars on eccen- 
tric, nearly Keplerian, approximately aligned orbits around a cen- 
tral black hole located at P2; PI is generated by orbital crowding 
of stars lingering at apocentre. Subsequent observations have been 
consistent with his model. Lauer et al. (1998) observed the nucleus 
with the Wide Field Planetary Camera 2 (WFPC2) on the corrected 
Hubble Space Telescope (HST), confirming the bimodal structure. 
Kormendy & Bender (1999) obtained spectroscopy with the Sub- 
arcsecond Imaging Spectrograph (SIS) finding an asymmetric rota- 
tion curve and a constant colour across the nucleus, showing PI and 
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P2's colours to be consistent with each other but not the bulge or 
a globular cluster. Further kinematics were recorded with the Faint 
Object Camera (Statler et al. 1999), the integral field spectrograph 
OASIS (Bacon et al. 2001) and the Space Telescope Spectroscopy 
Imaging Spectrograph (STIS) on HST (Bender et al. 2005). All of 
these observations show that the kinematic centre of the nucleus is 
very close to P2, as predicted by T95's eccentric disc model. 

Further refinement of this picture has come from studying the 
nucleus at ultraviolet wavelengths. Almost all of the UV emission 
from the nucleus comes from a tiny (< 0!'l) source located at P2 
(King et al. 1995; Lauer et al. 1998; Brown et al. 1998), whose 
optical-UV colours and spectra are consistent with a population of 
A stars (Lauer et al. 1998; Bender et al. 2005), distinct from the 
red K-type spectrum of the rest of the nucleus (Bender et al. 2005). 
Bender et al. (2005) find that this compact, young, blue population 
has a maximum velocity dispersion of a = 1183 ±201 km s~^, 
significantly higher than the red stars' ~ 250 km s~^. They label 
the UV peak P3 and have found that its photometry and kinemat- 
ics are well modelled by a separate, almost-circular disc of blue 
stars around a central black hole of mass (1.1-2.3) x 10^ M©. This 
provides very strong evidence in support of the presence of a su- 
permassive black hole, which is the most fundamental requirement 
of T95's model. Most recently, Lauer et al. (2012) has found that 
the surface brightness profile of the young population is described 
by an exponential profile of scale length Of'075 ±0f'01. 

Tremaine's original eccentric disc model consisted of three 
Keplerian orbits, coplanar with the main disc of M31, projected 
onto the sky and convolved with Gaussian point spread functions 
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(PSFs). The original model included neither disc self-gravity nor a 
realistic treatment of the disc's internal velocity dispersions but was 
still able to broadly reproduce kinematic features. To be stable over 
many dynamical times such a disc would require apsidal alignment 
to be maintained; T95 proposed that this could be achieved by the 
disc self-gravity and argued that two-body relaxation would lead 
to a disc thickness of 0.3 times the disc radius. Estimates of the 
disc mass from mass to light ratios place it at around IQp Mq. This 
is substantial enough to affect the dynamics of the disc although the 
nucleus falls within the sphere of influence of the black hole, which 
dominates the orbits and imposes regularity. A self-gravitating ec- 
centric disc mode maintaining orbital alignment would also pre- 
cess under a uniform pattern speed. Since the T95 model, several 
two-dimensional self-gravitating models with these properties have 
been constructed (Statler 1999; Bacon et al. 2001; Salow & Statler 
2001 ; Sambhus & Sridhar 2002; Salow & Statler 2004). These mas- 
sive models have found a variety of disc masses, pattern speeds and 
orbital distributions. 



Peiris & Tremaine (2003) took a different approach. They re- 
vised the 1995 model and constructed fully three-dimensional mod- 
els with a parametric distribution function that ignored the self 
gravity of the disc; the gravitational potential in their models is 
due solely to the central black hole, which greatly simplifies the 
modelling procedure. Their models are the most successful to date 
at fitting the observed kinematics. They also found that thickened 
disc models that were misaligned with respect to the large-scale 
M31 disc produced significantly better fits than coplanar models, 
echoing a result seen in the 2d models. 
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Figure 1. V-band image of the nucleus (Lauer et al. 1998) annotated to show 
the relative positions of PI, P2 and P3 in arcseconds. The pair of white lines 
represents the positioning and width of the STIS slit (Bender et al. 2005). 
The brightest point in the image is at PI and is marked with cross. Centred 
in the image and marked with a cross is the location of the black hole and 
P3. P2 is on the anti-Pl side of the black hole. 



The current picture of the nucleus (see figure 1) has the black 
hole (hereafter BH) embedded in P3, which is explained by a flat, 
circular exponential disc of blue stars Of' 033 from the photometric 
bulge centre. This young, blue disc is surrounded by the larger, 
red, eccentric disc. PI is made of stars crowding at apoapsis, while 
P2 is now identified as stars at pericenter in the elongated region 
on the anti-Pl side of P3. The system remains of great interest: 
the origin of the nuclei and the relationship between the red and 
blue populations are unexplained and the mass and pattern speed 
of precession of the disc have not been pinned down consistently. 
Understanding the dynamics of the disc will help determine the 
BH mass more accurately which is of use in better determining the 
relationship between a BH and the host galaxy. It is also of interest 
to confirm whether the eccentric disc is aligned with the main disc 
ofM31. 

In this paper we present a natural development of the mod- 
elling approach started by Peiris & Tremaine (2003). Like them, 
we ignore the self-gravity of the disc, leaving the construction of 
fully self-gravitating models for a subsequent paper. But instead 
of considering parameterised functional forms for the phase-space 
distribution function (hereafter DF) of the nucleus, we model the 
DF non-parametrically as a mixture of Gaussian rings whose am- 
plitudes are allowed to vary: one of our goals is to "let the data 
speak for themselves" and then examine closely the structure of the 
DF. The paper is organised as follows. In section 2 we summarise 
the data we use and additional post-processing applied. Section 3 
describes our modelling procedure and section 4 our results. Sec- 
tion 5 sums up. 

Throughout this paper we adopt a distance of 770kpc to M3 1 . 



2 OBSERVATIONAL DATA 

Our models use observations from three instruments: the HST pho- 
tometry recorded with WFPC2, detailed in Lauer et al 1998 (here- 
after L98), the high resolution kinematics from STIS (Bender et al 
2005, B05) and the kinematic maps from the integral field spectro- 
graph OASIS (Bacon et al 2001, hereafter BOl). 

2.1 Photometry 

We use the deconvolved V-band (F555W) WFPC2 photometry of 
L98. This image is 1024x1024 pixels in size with a pixel scale 
of Of'011375. P3 is located at (513, 518). We note that the orien- 
tation of the image has North 82.3° clockwise from the y-axis - 
the data has been rotated away from the original alignment of the 
CCD, so that the comers of the image are blank - but we make no 
attempt to reorient the image. This image has already been reduced 
and deconvolved but we perform some processing on the image in 
order to reduce this data to a tractable number of observable s for 
our fitting procedure. Bright foreground stars are first masked out 
and we crop the image toa512x512 pixel region centred on P3. 
This is a region of side 5f'824 and incorporates the nucleus and 
the inner parts of the bulge and removes the blank comers of the 
image. This cropped image is then rebinned to "super-pixels" of 
side 2^ pixels containing the mean of the 2^^ — n pixels, with n the 
number of masked pixels falling in the super-pixel. We consider 
two schemes for choosing super-pixels: in the first we simply take 
1 = 2 everywhere so that super-pixels contain a 4 x 4 pixel region. 
In the other / is chosen based on the surface brightness of the image 
scaled by r~^/^, where r is the radius from P3, to provide a good 
balance of detail around P2 and PI. This gives / = (the resolution 
of the deconvolved image) at PI and P2 and / = 4 in the outermost 
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parts of the image. A total of 4096 super-pixels are generated in 
this scheme. 

Ajhar et al (1997) and Lauer et al (1998) showed that that at 
the resolution of WFPC2 the underlying number of stars in M3 1 
shows strong surface brightness fluctuations. The limited number 
of stars is the dominant residual between the structure of the nu- 
cleus and any model we fit. We treat the noise from the fluctuations 
as the sole source of our errors in our photometry and ignore pho- 
ton noise so that the fractional error in a spatial bin is given by 
o = N~^^^, where N is the effective number of stars within the spa- 
tial bin. L98 gives the effective magnitude of each SBF "star" as 
fhi = 23.4, or 224 L©. The error associated with a spatial bin con- 
taining hLq is therefore V224nLQ. 

2.2 STIS kinematics 

The STIS kinematics are taken from Table 5 of Bender et al. (2005). 
These consist of the bulge- subtracted Gauss-Hermite coefficients 
derived with the Fourier Correlation Quotient Method at the 22 po- 
sitions listed for each of V, a, and h4 and their accompanying 
errors. We adopt the same slit widths of 0.1 ''and position angle of 
39°. The quoted positions are given with respect to the location of 
P3. 

We use equation (21) of Peiris & Tremaine (2003) to model 
the STIS PSF. When this is convolved with a a 0''1 -square top hat 
to represent the slit, the result is well approximated by a double 
Gaussian 



PSF(X, j) : 



li 



2naf 



(1) 



with amplitudes Ii = 0.24, 12 = 0.76 and dispersions Gi = 0f'042 
and 02 = Of'087 respectively. We assume that equation (1) as the 
effective PSF of the STIS observations. A more sophisticated treat- 
ment might take account of asymmetries in the PSF and the varia- 
tions in spatial binning along the slit. 

2.3 OASIS kinematics 

We also make use of the kinematics of the integral field spectro- 
graph OASIS, which were kindly provided by Eric Emsellem. We 
opt to use the higher S/N data set, M8. This consists of V , a, hi, and 
h4 values derived from spectra taken at 1123 positions, spaced by 
of' 09. We registered the image with the WFPC2 data using a sim- 
ilar process to BOl. We found a close match in angle (0.7°) and a 
small offset of (-Of' 02, -Of'03) between the two images. 

Like BOl we assumed that the OASIS measurements have a 
PSF that can be described by a sum of three Gaussians and allowed 
the parameters of these Gaussians to float freely in the registration 
process, taking the form 



PSF(X, j) : 



1 



3 



(2) 



The resulting PSF differs from that found by BOl with Gi = Of' 230, 
02 = Of'587, 03 = 0."440 and hjh = 0.836 and I3 /h = 0.057. 



3 MODELLING PROCEDURE 

Our models are straightforward generalisations of those of PT. In 
particular, we ignore the self gravity of the disc and the gravita- 
tional influence of the bulge and assume that the potential is purely 



Keplerian. The mass of the central black hole is the single free 
parameter in our model potential. We assume that the BH is located 
atP3. 



3.1 Coordinate systems 

Following PT03 we use three different coordinate systems: "orbit 
plane", "disc plane" and "sky plane". All three coordinate systems 
have origin O coincident with the BH. Our models have biaxial 
symmetry. This provides a natural definition of the "disc plane" 
{x,y,z) coordinate system: the model is symmetric under reflec- 
tions in the {x,y) and planes. The orbit-plane and sky-plane 
coordinate systems are defined as follows. In the potential of the 
BH all orbits are Keplerian ellipses. An orbit with semi-major axis 
a and eccentricity e defines a coordinate system (x',y,z') in which 
the Oz' axis is parallel to the orbit's angular momentum vector and 
the Ox' axis points towards pericentre. That is. 



X =a{cosE — e), 

y = a\/ l — e^ sinE, 
z' = 0, 



(3) 



where the eccentric anomaly E is related to the mean anomaly M 
through Kepler's equation M = E — esinE. In this coordinate sys- 
tem the apocentre of the orbit is located at (x',y,z') = (— a(l + 
^),0,0) and the pericentre at (a(l — ^),0,0). Each star defines its 
own orbit-plane (x',y,z') coordinate system. 

A star's disc-plane coordinates {x,y,z) are related to its orbit- 
plane coordinates (x',y,z') through 




cos CI —sinQ. 
sin 12 cosQ. 
1 

(cos CO —sin CO 
sino) cos CO 





(4) 



where the angles ca, / and Q. are the star's argument of periapsis, the 
inclination and longitude of the ascending node, respectively. Any 
orbit in the Keplerian potential of the BH can be labelled by the 
five integrals of motion (a,^,ca,/,^2), but it proves convenient to 
replace e and co by the eccentricity vector e = (^cosG5, ^sinG5, 0) = 
{ex,ey,0), where the longitude of periapsis 05 = CO + 11. The vector 
e points from the BH towards the projection of the pericentre onto 
the z = disc plane. Its magnitude is the scalar eccentricity e of the 
orbit. 

Projected, sky-plane coordinates (X,7,Z) are related to disc 
plane coordinates via 

— sin 9/ 
cos 9/ 






cosQa — sin9<3 
sin 9^ cosQa 





(5) 



in which (9^,9/, 9/) are the three Euler angles specifying the ori- 
entation of the disc with respect to the observer's reference frame. 
The (X, Y) plane is the sky plane, with the positive X axis pointing 
west and positive Y axis north. The Z axis then points towards the 
observer; the line-of- sight velocity is therefore Vios = — Z, follow- 
ing the usual convention that that receding objects have Vios > 0. 
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3.2 Distribution function 

We model only the old stellar population of the disc; the young 
stars around P3 and the bulge are treated as contaminants (see §3.4 
below). Assuming that the old stars are collisionless and homoge- 
nous, their dynamics can be completely described by a distribution 
function (DF) /(x,v), defined such that /(x,v)d^xd^v is the ex- 
pected number of such stars within a small phase- space volume 
d^xd^v around the point (x,v). We further assume that the system 
is in a steady state. By Jeans' theorem the DF must be expressible 
as a function of the integrals of motion (a, e, 7,11). The number of 
stars within a small phase-space volume element d^xd^v is then 

/(x,v)d^xd\ = /(a,e,/) ^(GM.)^/^^^/^ sin/d/dadedMM, (6) 

the \{GM.fl^ a I mil factor comins from the Jacobian relating 
(x,v) to (a,e,/,Il,M). 

We assume that the distribution function can be decomposed 
into a weighted sum of rings, 



j 



(7) 



in which each ring has a uniform distribution of II G [0, 2n) and 
Gaussian distributions in a, e, /, 



fj(a,e,I) =NjQxp 



(a — aj)'^ 






[ 2a,,/ J 


exp 


[ 2ae,/ J 



(8) 



exp 



2<^/,/ 



and it is understood that fj{a,e,I) = if a < or |e| > 1. The nor- 
malisation factor Nj is included to give each ring unit total lumi- 
nosity. We use a Monte Carlo method to construct the rings, which 
avoids explicit calculation of Nj. 

Our rings span a (20 x 9 x 4) grid in {aj,exj,Oij), with 
mean semimajor axis aj running logarithmically from to 
10'^ mean x component of eccentricity vector exj from —0.8 to 
+0.8 in steps Ae = 0.2, and the dispersion in inclination ajj 
drawn from {12°, 24°, 36°, 48°}. For the spreads in {a,e) we take 
<^aj = 0.8ayAloga and Oej = O.SAe = 0.16. All rings have mean 
Cyj = 0. Together with the zero mean inclination of each ring this 
makes our models symmetric under reflection in the j = and z = 
planes. Rings with mean ex>0 are "aligned" in the sense that their 
apocentre (and therefore peak density) lies somewhere along the 
negative x axis. Rings with mean ex <0 have the opposite orienta- 
tion. 



3.3 Observables of each ring 

We compute the projected properties of each ring fj (equ. 8) 
by sampling it with 10^ points, drawing A^^ = 5000 values of 
(a,e,/,Il,co) from the distribution /y(a,e,/)a^/^ sin/ that appears 
on the right-hand side of (6), then sampling 200 points equispaced 
in mean anomaly M along each of these orbits. To avoid sampling 
artefacts we select the initial value of M for each orbit at random. 



3.3.1 WFPC 

We use this sample of 10^ points to estimate the projected surface 
density (or, equivalently, the zeroth velocity moment) 



: J dZdVxdVydVzfj 



(9) 



on a fine grid of Of'011375 x Of'011375 pixels on the (X,7) sky 
plane by simply counting the number of points that fall within each 
pixel. The contribution of the ring to each WFPC2 "superpixel" 
is found by taking the mean value of A/7""''(X,F) over the relevant 
range of {X,Y). We do not carry out any PSF convolution for the 
WFPC2 photometry because the image of Lauer et al. (1998) has 
already been deconvolved. 



3.3.2 OASIS 

The OASIS kinematics of Bacon et al. (2001) consist of measure- 
ments of mean velocity V{X,Y) and velocity dispersion a(X,7), 
together with Gauss-Hermite coefficients /z3(X,7) and h4{X,Y). 
We ignore their and h4 measurements and assume that their 
measured V(X,Y) and (V^ + a^)(X,7) distributions probe directly 
the (PSF-convolved) first- and second-order moments of the line- 
of-sight velocity distribution (LOSVD) at the point (X,7) on the 
sky. This is a good approximation provided the underlying PSF- 
convolved LOSVDs are reasonably close to Gaussian. For mod- 
elling purposes it is more natural to consider luminosity- weighted 
moments 



IjHx,Y)=I{X,Y)V{X,Y), 
l/{X,Y)=I{X,Y){V^ + G^){X,Y), 



(10) 



where I(X,Y) is the underlying surface brightness. We ob- 
tain I{XJ) by convolving the WFPC2 image with the OASIS 
PSF (2). The observational errors on these first- and second- order 
luminosity- weighted velocity moments are obtained by adding the 
uncertainties on /, V and a in quadrature in the obvious way: 



(Aa/ 



1\2 . 



{AI)^V^+I^{AVY 



{Ai?f = (A/)2(y2 + 0^) +f{2\V\AVf +/2(2aAa)^ 



(11) 



Following (9) above, the contribution of the ring to the 
luminosity- weighted first and second moments of the line-of- sight 
velocity distribution. 



A/].'^^'(X,7): 



2,fine 



: I dZdVxdVYdVzVzfj, 
: I dZdVxdVydVzVifj, 



(12) 



are estimated by weighting each of the 10 sample points by Vz and 
yj, respectively. Having these jLij^^^{X,Y) distributions we con- 
volve with the OASIS PSF (2) to obtain the contribution the ring 
makes to the model's predictions for the first and second velocity 
moments of the OASIS kinematics. 



3.3.3 STIS 

Our first attempt at fitting the STIS kinematics was based on the 
same assumption that we could use the STIS V and a measured 
by B-i-05 as direct estimates of the (STIS PSF-convolved) first- and 
second-order moments of the LOSVD. That did not work well: at 
STIS resolution the LOSVDs are far from Gaussian, as we show 
on Figure 2, and the strong high- velocity wings caused by the BH 
mean that it is not possible to obtain reliable estimates of the first 
and second moments from the observed spectra. This last point was 
one of the motivations for the introduction of Gauss-Hermite se- 
ries (van der Marel & Franx 1993) to parametrize LOSVDs. There- 
fore we fit our models to B-i-05's Gauss-Hermite parametrizations 
of the LOSVDs along the STIS slit. Our method for fitting the 
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Figure 2. LOSVD of a typical model at offset R = 0'.'2 along the STIS slit 
(solid blue curve) and its 4-th order Gauss-Hermite reconstruction (equ 16, 
dashed red curve). The best-fit Gaussian to the LOSVD has (y,y,o) = 
(0.81, — 85kms~\281kms~^). For comparison, the classical zeroth-, first- 
and second-order moments of this LOSVD are 1, Vq = — 36kms"^ and 
Go = 759kms~\ the latter being dominated by the strong high-velocity 
wings. 



Gauss-Hermite coefficients follows the same lines used in other 
orbit- superposition models (e.g., Cretton et al. (1999)), but taking 
extra care to treat the normalisation of the LOSVD s correctly. 

We recall some details of Gauss-Hermite expansions. Sup- 
pose that we are given an LOSVD Lq{v), normalised such that 
/ Lo(v) dv = L The Gauss-Hermite expansion of this Lo(v) is 



^ 7=0 



(13) 



where w = (v — V)/a, a(w) = / \^2n is the standard Gaus- 

sian and the Hj(w) are Hermite polynomials. We adopt vdMF93's 
normalisation for the latter. Using the orthogonality properties of 
the Hj, it is easy to show that the Gauss-Hermite coefficients hi are 
given by 



/^/(y,V,a): 



y 



/: 



L{)(y)Hi(w)a(w)dv. 



(14) 



That is, there is a different Gauss-Hermite series (13) for each 
choice of (y, V,a). The "line- strength" parameter y simply scales 
all the hj, but it proves important as we shall now see. 

A particularly natural choice of the parameters (y,y,a) are 
those that minimise 

m2 



Xo ' 



/: 



io(v) 



Ya(H') 



dv, 



(15) 



in which case it can be shown that the first few Gauss-Hermite 
coefficients from (14) become = \ and /^i = = 0. In other 
words, if we choose (y,y,a) to be the parameters of the best-fit 
Gaussian to the LOSVD L(v), then the LOSVD can be written as 



L(V): 



ya(w) 



7=3 



(16) 



with /z3, /z4, ... given by the integral (14). Conversely, if we adopt 
the parametrization (16) and fit (y, V, a, /^s , /i4 , . . . ) simultaneously to 
Lo(v), then we get back the same parameters we would obtain by 



first fitting (y,y,a) by minimising (15) and then using (14) to find 
the hi. For a strongly non-Gaussian Lo(v) the parameters (y,y,a) 
obtained by minimising (15) need not be close to zeroth-, first- and 
second-order moments of Lo(v), as shown on Figure 2. In particu- 
lar, y need not be close to 1 . 

Gauss-Hermite fits to the LOSVDs of real galaxies, includ- 
ing B-i-05's measurements of M31, adopt the parametrization (16) 
and fit (y,y,a,/i3,/i4, ...). Unfortunately, the parameter yis rarely 
reported, presumably because it is strongly affected by systematic 
effects in the fitting procedure, such as template mismatch, and be- 
cause it does not affect the shape or width of the LOSVDs. Never- 
theless, we note that it is an essential part of any Gauss-Hermite ex- 
pansion. For now we assume that y is known. Our iterative scheme 
for reconstructing it from the models is described in section 4.3 
below. 

Notice that equation (14) shows that the Gauss-Hermite co- 
efficients hi can be thought of as modified moments: each hi is the 
integral over velocity space of some linear combination of the clas- 



sical moments l,w,w 



but weighted by the Gaussian fac- 



tor a(w). Therefore, given a Gauss-Hermite fit (y,y,a,/i3,/z4) to 
the line-of-sight velocity distribution at projected position (X,F), 
we treat (y, V, a) as being known perfectly and take the (luminosity- 
weighted) modified moments (A^obs' •••'^obs) = {^i^i^ih^^^h^) 
as our observables, where the surface brightness /(X,F) is ob- 
tained by convolving the WFPC2 image by the STIS PSF (1). 
We use equations (10) of vdMF93 to propagate B-i-05's quoted 
uncertainties (AV/, Ag/, A/z3 A/z4 /) to our observational errors 



'Mobs) 



Just as for the classical velocity moments, for each ring j = 
1,2, ... we use our Monte Carlo sample of 10^ positions and veloc- 
ities to evaluate the expressions 

f/ 

1 



Aif"^(X,F|Y,V,-,0,-; 



dZdVxdVydVz 



X exp 



(17) 



Hk(w)fj 



for the modified moments, in which the rescaled velocity 

_ {-Vz)-Vi 



(18) 



Then we convolve each of these distributions with the STIS PSF (1) 
and read off the values at = giving the contribution 

of the 7* ring to the modified moment of the STIS data point, 

^?,obs' 

3.4 Modelling the effects of the bulge and P3 

Our ring system is designed to model only the old red stars of the 
eccentric disc, but some of the light observed in the central few arc- 
sec of M31 comes from other sources. The two main contaminants 
are M3rs bulge and the compact young stellar cluster at P3. 

We follow Kormendy & Bender (1999) in modelling the 
surface brightness of the bulge as a Sersic profile I{R) = 
/oexp(-(i?//?„)i/") ^-^j^ -j^^g^ ^ ^ 2.19, scale radius Rn = 14f'0 
and central V-band surface brightness Iq = 15.40 mag. Our kine- 
matic model for the bulge is very simple: it is non rotating and 
has a constant velocity dispersion Gbuige = 120km s~ ^ . We add the 
contribution from this simple bulge model to our models' predic- 
tions for the WFPC2 photometry and the OASIS V and a maps. 
We do not add it to the STIS predictions; we assume that B05 have 
successfully removed the bulge contribution from their STIS kine- 
matics. 
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To model the contribution the young stars from P3 make to 
the V-band light we include a another component having surface 
brightness I.{X,Y) = E3 exp(-/?//?o) in which =X^^Y^ and 
the scale length Rq = 0'/015 (L12). We follow B05 in assuming that 
the kinematics extracted from the red spectra are unaffected by the 
young stars; these stars affect only the WFPC2 photometry, not the 
OASIS or STIS kinematics. 



3.5 Fitting the weights 

For given BH mass and disc orientation (9<2, 6/, 9/), the model's 
prediction for any observable 0/ can be be written as Y^jPij^j^ 
in which wj is the weight given to the 7* ring (equ. 7) and the 
matrix Pij gives the contribution that the 7* ring makes to the 
observable, calculated using the method described in §3.3 above. 
An observable 0/ can be the light within a "superpixel" (given 
by WFPC2 photometry), a classical first- or second-order velocity 
moment (OASIS kinematics) or a 0*...4*-order modified moment 
(STIS kinematics). We do not include the zeroth-order moments of 
the OASIS kinematics as these contain no additional information 
over the WFPC2 photometry. 

Having a vector of observables (Oi,02,....) and associ- 
ated uncertainties (Ai, A2, . . .), we use a non-negative linear least 
squares algorithm (Lawson & Hanson 1974) to find the vector of 
non-negative weights w that minimises 

.2_v {Oi-'LjPij^P^ 



where = 0.307. The parametric form for the backbone eccen- 
tricities is given by 



A/ 



(19) 



We then take this minimum value of as a measure of the good- 
ness of fit of the model with parameters (M., 9/, 9/, 9^2). 



4 RESULTS 

Ideally, we would like to include all data sets (WFPC2, STIS, OA- 
SIS) in our set of 0/ and fit simultaneously to them. However, cal- 
culating the contribution the 7* ring makes to a particular modified 
moment is computationally intensive and a full study of the param- 
eter space spanned by M., 9/, 9/ and 9^ that includes the STIS and 
OASIS data sets is not viable. We therefore conduct our modelling 
in two stages. We first fit models to the WFPC2 photometry and 
the OASIS velocity distribution in order to find the best set of ori- 
entation angles (9^,9/, 9/). Then, having these angles, we fit to the 
WFPC2 photometry and STIS kinematics to obtain our estimates of 
the BH mass and the structure of the phase- space DF of the nucleus. 
As a final test of this model, we "observe" it at OASIS resolution 
and compare it (by eye) to the real M3 1 . 

Before embarking on any of this model fitting, however, we 
first confirm that our models are indeed able to reproduce some of 
PT03's results. 



4.1 A test: reproducing PT03's model 

An immediate test of our modelling machinery is to reproduce 
some of PT03's results by using a sum of rings (8) that approx- 
imates one of their DFs. We focus here on their favoured non- 
aligned model, but we found comparable results with the poorer- 
fitting model that is forced to be aligned with the main M31 disc. 
PT03's non-aligned models have DF 

[e-e;^(a)]^ ' 



em(^) = oc(ae — a) exp 



{a-agf 



2w2 



(21) 



/(a,e,/) =g(a)exp 



2a2 





r fi 1 


jexp 


2a/(a)2_ 



(20) 



with a = 0. 197 pc"^ = 4.45 pc, ag = 1.71 pc and w= 1.52 pc. 
For the dispersion in inclination the form is 

a/(a) = exp(-a/a/) (22) 

where = 24? 6 and a/ = 31.5 pc. The function g(a) that sets the 
semi-major axis distribution is 

/ X y (2^exp(-a/ao) 

^^""^ 'l + exp[ci(a-C2)] 27r2(GM.)3/2a2a2(a)' ^ ^ 

withao = 1.37 pc,M, = 10.2 x lO^M©, ci =4pc"^ and C2 =4.24 
pc. We treat the overall normalisation Zq as a free parameter. 

We approximate this DF by a sum of 720 rings of the form (8) 
whose semi-major axes aj are spaced logarithmically in radius 
from 0.1 pc to 19 pc, with Oaj = O.SajAloga. We set a^j = 0.307 
for all rings and set the mean eccentricity ej and dispersion in in- 
clination Gij of the 7* ring by evaluating (21) and (22) above at 
a = aj. The weights wj are set proportional to ajg{aj), the factor 
of aj coming from the fact that our rings are equispaced in log a. 

Figure 3 shows the resulting model predictions for the WFPC2 
photometry when viewed at the same angles as PT03's best-fit non- 
aligned model = -34?5, 9, = 54? 1, 9/ = -42.8P2) and includ- 
ing the contribution of the bulge model from sec. 3.4. Our recon- 
struction of the projected surface brightness in their model agrees 
closely with their figure 3: the model produces a nucleus that is 
broader than the observations and has an overly extended flat pro- 
file at PI. 

4.2 Determining the orientation of the disc 

We conduct an exhaustive scan over the space of orientation angles 
(9(3, 9/, 9/) and black hole mass to obtain our best guess for the 
orientation of the disc. We first do this for a model that is forced 
to be aligned (9/ = 77?5) with the larger-scale M31 disc, before 
letting 9/ vary freely. We include in these fits the 1123 OASIS V 
data points as well as a broad field spanning 5^^824 x 5^^824 from 
the WFPC2 data: this utilises our varied super-pixel scheme and 
contains 4096 data points. We do not include the OASIS a maps in 
the fits, as it is dangerous to assume that the measured + is a 
good estimate of the true second moment of the LOSVD (figure 2). 

For the aligned model we fix 9/ = 11? 5 and take and 9/ in 
the intervals [-90°, 0°] and [-55°, -25°] with spacing A9 = 2°. 
A contour plot depicting the shape of the goodness of fit for 
the inner region of this range is shown in figure 4. As expected the 
aligned model constrains 9/ tightly but there is a weak dependence 
on 9(3, with the best fitting model falling at = —60° and 9/ = 
-35°. 

The non-aligned model looks at 9^ from —58° to —20°, 9/ 
from 50° to 69° and 9/ from -39° to -29°, with A9« = 2° and 
A9/ = AQi = 1°. The best fit appears at 9« = -34°, 9/ = 57°, 
9/ = -35° for a black hole mass M, 1.25 x lO^M©. The shape 
of the goodness of fit for the non-aligned model is shown in 
figure 5. The presented scan is for a single random distribution 
of stars projected with different angles: scans using distributions 
drawn from a different random seed found broadly similar results 
with slight variations (figure 6). This presents issues with deter- 
mining a single best fitting set of angles. Using a much larger (10^) 



© 2013 RAS, MNRAS 000, 1-13 



Nucleus ofM31 1 




-1 1 

X (arcseconds) 




-1 

X (arcseconds) 



Figure 3. Nuclear V-band surface brightness distribution. The left panel shows the data and the right our reconstruction of Peiris & Tremaine's (2003) 
non-aligned model. Contours are at 0.25 mag intervals. Compare to figure 3 of PT03. 




Figure 4. for 0/ vs. for the aligned model. Contours are spaced at 
Ax^ = 4 intervals. 



number of stars per disc in the region around the best fit angles in- 
formed our final selection of specific angles, however even then the 
random seed affected results. While in principle the angles could be 
determined more accurately, in practice shot noise from the model 
(and also the finite number of stars in the real nucleus!) limits 
the accuracy in each angle to the order of 1°. Our final selec- 
tion of angles was determined from the projection of the likelihood 
exp(— ^X^) onto each of the 9^, 9/ and 9/ axes. 

We have experimented with allowing the bulge surface bright- 
ness /q and the central surface brightness E3 of the young stellar 
disc described in sec. 3.4 to float our fitting procedure by including 



ditional columns to the projection matrix P/y, but we find that this 
makes little difference to our results. 



4.3 Fitting WFPC2 photometry and STIS kinematics 

Having the orientation angles we now drop the OASIS V maps 
and focus on using WFPC2 photometry together with the STIS 
LOSVDs to further constrain the model. The region of the WFPC2 
photometry we use is restricted to an ellipse of semi-major axis \ '.'6 
and axis ratio 0.6 centred on P2. This ellipse is just large enough 
to encompass all STIS positions. We rebin Lauer et al. (1998)'s 
dithered image 4 by 4 into "superpixels" of side Of' 0455. Our vec- 
tor of observables 0/ consists of the WFPC2 fluxes in all 2335 such 
superpixels that lie within the ellipse, together with the 5 modified 
moments (^/^o 5 • • • 5 ^/,4 ) obtained from (y/ , V/ , a/ , /z3 , h/^. ) for each of 
the / = 1...22 LOSVDs measured by Bender et al. (2005) using the 
procedure described earlier in section 3.3.3. 



4.3.1 Reconstruction of the y{R) profile 



them as additional "weights" wj in the model and adding two ad- 



The only complication in this is that we do not know the line 
strengths y for any of our LOSVDs. We do, however, expect y to 
be reasonably close to one and so we first fit a model in which all 
ji = 1. Then, having the weights wj, we construct a realization of 
this model and "observe" it convolved with the STIS PSF (1). We fit 
Gauss-Hermite coefficients (y,y,a,/z3,/i4) to the model LOSVDs 
at each point along the STIS slit, giving us a more informed esti- 
mate of how y varies along the slit. We then repeat the whole fitting 
procedure, replacing our original y/ = 1 guesses with values read off 
from this reconstructed y{R) distribution. We find that the resulting 
model converges after only a couple of iterations of this scheme. 
Figure 7 plots representative y profiles obtained by this procedure 
for models with black hole masses M^/IO^M© = 1.0, 0.8 and 1.2. 
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Oi 

Figure 5. Contour maps of for the slice through the best fitting model. From left to right: 0/ vs. 0^ for 0/ = -35°; 0/ vs. 0^ for 0/ = 57°; 0/ vs. 0/ for 
= -34°. Contours are spaced by A^^ = 25. 




Figure 6. Zoomed in contour maps. From left to right: 0, vs. 0^^; 0/ vs. 0^; 0/ vs. 0/. Contours are spaced by Ax^ = 7.5. The black contours show the same 
distribution as figure 5. The grey contours show an alternative realisation of the distribution function to illustrate the effect of shot noise. 
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Figure 7. The Gauss-Hermite "line- strength" parameter y reconstructed 
along the STIS slit using the iterative procedure described in section 4.3 
for three different assumed black hole masses, as indicated on the legend. 



This shows that y is significantly depressed close to the black hole 
where the LOSVDs are least well described by simple Gaussians. 



4.3.2 Comparison of best-fit model against observations 

The formal of the model with M. = 1.0 x lO^M© is 782. To 
put this in perspective, the vector of observables 0/ that the model 
fits has 2445 elements: 2335 WFPC2 fluxes plus 5 x 22 modified 



moments. For comparison, the model with = 0.8 x 10^ has 
= 825, while the model with M, = 1.2 x 10^ M© has = 859. 
Figure 8 shows the Gauss-Hermite coefficients of the reconstructed 
models along the STIS slit. The agreement with B05's observed 
kinematics is good, but there some features the best = 1.0 x 
10^ Mq model cannot reproduce: the model does not fit the detailed 
shape of the V{R) and h^iR) profiles between R = -0'.'9 and Of'3 
well and it predicts a central a{R) that is slightly too high. 

Figure 9 shows how well the = 1.0 x lO^M© model fits 
to the WFPC2 photometry. Results for the other two black hole 
masses are similar. The agreement is good in the central regions, 
but beyond about 1 arcsec from P2 the model's surface brightness 
profile falls off too slowly compared to the observations. One pos- 
sible explanation for this is that our model is simply too coarse; 
the thickness of each of the na = 20 rings (equ. 8) in these models 
is a a = 0.23a, which sets the models' characteristic radial spatial 
resolution. Another is that our model for the contribution of the 
bulge light (sec. 3.4) may be wrong within the innermost couple of 
arcsec. 

Based on these comparisons, we interpret the relatively low 
values of x^ of our models not as a indication of the outstanding 
quality of our model fits, but instead as a sign that the treatment of 
the observational uncertainties - particularly of the WFPC2 pho- 
tometry - could be improved. Nevertheless, we believe that the 
reader will agree that simple "chi-by-eye" tests indicate that our 
models produce the best fits to date of the M3 1 eccentric disc sys- 
tem. Of course, this is to be expected given that we have > 700 free 
parameters to play with, which is at least an order of magnitude 
more than most previous models of M3rs nucleus. 
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Figure 8. Gauss-Hemiite coefficients parametrising the LOSVDs along the STIS sHt measured by B05 (points) together with our model fits for three different 
assumed black hole masses (curves). 
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Figure 9. Observed V-band WFPC2 image (blue) with the fit from our 
M. = 1.0 X 10^ Mq overlaid on top (red). Contours are spaced at 0.3 mag- 
nitude intervals. 

4. 3. 3 DF of best-fitting model 

What can we learn from all these free parameters, specifically 
the orbit weights wp. Figure 10 shows two views of the DF 
© 2013 RAS, MNRAS 000, 1-13 



f{a,ex,eyj) of our best-fit model. This model has a strong neg- 
ative eccentricity gradient between 0'/5 and l'/2 (corresponding to 
PI), which is very similar to PTOS's best-fit model; even the dis- 
persion in eccentricity (figure 1 1) is similar to their value of 0.307. 
There are three important differences between our model and theirs, 
however. 

(i) Beyond the mean eccentricity in our models becomes 
negative, meaning that the rings become mildly antialigned. We 
find that the strength of this feature depends on the details of our 
bulge model and so it is hard to judge its significance, but we note 
that just such a feature was predicted by Statler (1999) in his anal- 
ysis of thin, self-gravitating discs. 

(ii) Whereas PT03's parametrised model had an exponentially 
declining a/ profile, our model fits a a/ profile that increases 
with radius for a > Of' 15. This last point is qualitatively consis- 
tent with the predictions of collisional models of disc evolution 
(Stewart & Ida 2000; Peiris & Tremaine 2003). The detailed agree- 
ment is not so good though: whereas the collisional models predict 
Gi/Ce — 0.5, our models fit Ci/Ce — 2. It is not immediately clear, 
however, how far one can apply these calculations that assume al- 
most circular e = discs to the strongly eccentric disc in M3 1 . 

(iii) f{a) in our model falls steeply towards the centre from a = 



10 C. K. Brown and J. Magorrian 




Figure 10. The distribution of orbits in our best-fit model with M. = 1.0 x 10^ M©. Left panel: the projected DF f{a,ex) obtained from the full DF f{a,ex,ey,I) 
by averaging over ey and /. The yellow curve plots the a(ex) profile of PT03's best-fit model for comparison. Right panel: RMS inclination angle (in degrees) 
as a function of {a,ex) obtained from /(a^ex^ey,!) by averaging over Cy. The f(a,ex) distribution from the left panel is overlaid as contours. 
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Figure 11. Mean eccentricity ex(a) and dispersions in eccentricity Ce^x{<^) 
and inclination Gi{a) as a function of semi-major axis a. Unlike figure 10 
the dispersion a/ in this plot is given in radians, not degrees. 

0^/5 to a ~ Of' 15, inside which there is antialigned (ex — —0.5), 
fat (a/ = 48°) distribution of orbits. Recall that the range of a/ 
reproducible by our chosen sample of rings is from 12° to 48°. 
This distinctive change in the distribution of the red stars is almost 
cospatial with the young. A- star population that make up P3. 

4.3.4 Other projections of the best-fit model 

Figure 12 shows the face-on and edge-on projected density distribu- 
tions of our best-fit model. Our machinery fits a density distribution 
which, when viewed face on, is broadly similar to the density dis- 
tribution adopted by Salow & Statler (2004) in their self-gravitating 
razor- thin models. Apart from our neglect of self gravity, the most 
significant difference between our models and theirs is that ours 
have a secondary density peak around P2 and also have significant 
thickness. 

Although the models we present in this section are not fit to 
the OASIS maps, we can nevertheless compare our reconstructed 



models' predictions against the real OASIS maps. Figure 13 shows 
the results for the = 1 .0 x 10^ Mq model. The model reproduces 
the shape and orientation of velocity and dispersion maps very well. 
This provides an independent test of the orientation angles and the 
broad-brush features of the DF inferred by our models. 

4.3.5 A brief experiment: retrograde orbits 

Finally, we note that one proposed origin of the eccentric disc in 
M3 1 is from an m = 1 instability in an initially circular stellar disc 
that contains a small counter-rotating population of stars (Touma 
2002; Jog & Combes 2009; Kazandjian & Touma 2012). A^-body 
simulations of this instability (Kazandjian & Touma 2012) show 
that the dominant, prograde population becomes the eccentric disc, 
while the minor, retrograde population is puffed up into a strongly 
triaxial distribution. As a quick experiment to test whether this is 
easily detectable, we have tried doubling up our ring distribution 
by including in our models the retrograde ring corresponding to 
each of the 720 prograde rings considered above, giving a total 
of 1440 rings. This produces a noticeably better fit to the obser- 
vations, with = 784 instead of 812, with the fitting procedure 
picking out three main counter-rotating components (figure 14) at 
{a,e) = (Of'2,0.8), (lf'2,0.6) and (3,-0.1). None of these is easily 
identifiable with the more diffuse component predicted by Kazand- 
jian & Touma (2012), although we suspect that that might be hard 
to detect with naive models such as ours. We interpret the inner- 
most retrograde ring as suggesting that our ring-based decomposi- 
tion fails at radii r < Of' 2, perhaps by not being fine enough. Simi- 
larly, we take the presence of the other two rings as a hint that our 
bulge model could perhaps be improved. 



5 CONCLUSIONS 

We have constructed an eccentric disc models of the nucleus of 
M31 by modelling it as a linear combination of fattened rings of 
stars orbiting in a purely Keplerian potential. Our models are an 
obvious generalisation of those of PT03 and - as expected - are 
noticeably (albeit perhaps not significantly) more successful at re- 
producing the features observed in the inner arcsec or so of M31. 



© 2013 RAS, MNRAS 000, 1-13 



Nucleus ofMSl 1 1 





Figure 12. Face-on (xy, left) and edge-on (xz, right) projected surface brightness distributions of our best-fit model. Contours are spaced at 0.2 magnitudes. 





Figure 13. OASIS V (left) and o (right) maps predicted by our models (white contours) compared to the observations from BOl. The contours in both panels 
are spaced 25kms~^ apart. The V contours run from to ±175kms~\ a from 150 to 225kms~^. 



One of the most fundamental parameters of any model of this 
system is its orientation. Like PT03, we assume that the disc has 
biaxial symmetry, but we find some differences in the Euler angles 
(9«, 9/, 9/) that specify the orientation of these symmetry axes. Our 
values of 9^ = —34° and 9/ = 57° agree reasonably well with their 
da = -34?5 and 9/ = 54.°51, but our value of 9/ = -35° (which 
directly controls the position angle on the sky) differs significantly 
from their 9/ = —42.° 8, more than can be accounted for by the ef- 
fects of shot noise in either the models or in the distribution of stars 
in the real disc. 

Our orientation is also consistent with the N-body model of 
BOl (9/ = 55° ± 5° and 9/ = -36°) and the models of Salow & 
Statler (2004) (9, = 63?51 ± 10?80 and assumed 9/ = -33?6) and 
Sambhus & Sridhar (2002) (9, = 51P54 and 9/ = 27.°34), though it 
should be noted all these models are 2d which imposes constraints 
on their geometry; Sambhus & Sridhar (2002) obtained their ori- 
entation by de-projecting the photometry of the disc such that the 
outer isophotes became circular. We note that these values of 9/ 



found for the old, red distribution of stars are close to the inclina- 
tion of the young population in P3 (9/ = 55° ±2°) measured by 
B+05. 

The most interesting result of our models is the DF they infer 
from the data, bearing in mind that they contain no prior "wisdom" 
about which DFs are dynamically plausible. The models suggest 
the presence of a distinct, compact disc of red stars within Of' 15 
of the black hole. This disc is anti-aligned (with eccentricity e < 
0) with respect to the larger-scale eccentric disc (which has e > 
0). Outside this compact region the eccentricity distribution of our 
models is very similar to PT03's. The main difference between our 
models and theirs is in the inclination distribution: our models fit 
an RMS inclination profile a/ (a) that closely tracks the dispersion 
in eccentricity, with ai{a)/ae{a) ^ 2. This is interesting: models 
of the collisional evolution of (circular) discs predict that this ratio 
should be ^ 0.5. 

Another difference between our models and PT03 is that we 
find evidence for an anti-aligned feature < 0) at a > I'/ 5. How- 
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Figure 14. The distribution of orbits in a model that includes both prograde and retrograde orbits. The top two panels show the DF f{a,ex 
inclination (right) for the prograde stars. The bottom two panels show the corresponding distributions for the retrograde population. 
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ever, the details of this feature depend on our assumed bulge model, 
which merits further investigation. 

Our models prefer black hole masses of the order of 1 .0 x 
lO^M©; masses higher than about 1.2 x 10^ are weakly ruled 
out. Although it would be possible to use our machinery to carry out 
a full scan of black hole masses and orientation angles, followed by 
a systematic investigation of the degeneracies in the DF, we believe 
that a more pressing task is to include the self gravity of the disc. 
The stars contribute a significant fraction (~ 20%) of the mass of 
the BH+eccentric disc system, which means that it is dangerous 
to read too much into our present, purely Keplerian models. Past 
2d (Bacon et al. 2001; Sambhus & Sridhar 2002; Salow & Statler 
2004, e.g.,) models have shown it is possible to get plausibly good 
fits to the kinematics for large disc masses. The space of 3d disc 
distributions that project to yield the observed surface brightness 
profile is degenerate, but our orbital ring system serves as a good 
starting point for this investigation. Work on self-gravitating self- 
consistent 3d disc models is now underway, and we hope that such 
models will provide further insight into the origin of the eccentric 
disc in M31 and possibly elsewhere (Lauer et al. 2005). 
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